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ABSTRACT 


Numerical  solutions  for  laminar  incompressible 
fluid  flows  past  thin  elliptic  cylinders  at  various  angles 
of  attack  are  obtained  using  a  stream  function-vorticity 
formulation.  The  flow  is  considered  to  be  two- 
dimensional  and  time-dependent.  Potential  flow  is 
selected  as  the  initial  condition.  Almost  steady-state 
solutions  have  been  computed  for  the  cases  a=  0°, 

Re  =  200;  45°,  Re  *  15;  and  a  =  90°,  Re  =  10 

where  a  and  Re  are  the  angle  of  attack  and  the 
Reynolds  number,,  respectively.  Vortex  shedding  has 
been  studied  for  a  =  45°s  Re  ~  30  and  200. 
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INTRODUCTION 


Flows  past  a  finite  flat  plate  of  infinite  span  in  a  channel  or  in 

an  infinite  body  of  fluid  have  been  studied  numerically  by  a  number  of 
1  -7 

investigators  „  In  all  these  studies  the  angle  of  attack  ex  is  either 

q  q 

0  or  90  except  for  work  reported  in  a  recent  paper  ,  in  which  the 

steady  flow  past  an  infinitely  thin  plate  between  two  parallel  walls 

o  2 

with  a-  45  was  investigated  Only  the  work  by  Fromm  and  Harlow 

considers  vortex  shedding  from  a  plate.  The  angle  of  attack  is  a  -  90° 

in  this  case. 

Periodic  vortex  shedding  from  a  circular  cylinder  has  been  the 

8-11 

subject  of  recent  numerical  studies  „  Extensive  surveys  on 

experimental  as  well  as  theoretical  work  on  periodic  vortex  shedding 

12  13 

from  bodies  were  presented  by  von  Krzywoblocki  and  Morkovin  . 

Numerical  studies  on  flows  around  plates ,  which  have  high  surface 

curvature  at  the  tips,  pose  many  more  difficulties  and  demand  much 

more  computer  time  than  studies  on  motions  past  circular  cylinders. 

The  analytical  and  numerical  problems  involved  have  been  investigated 

14-17 

and  discussed  elsewhere  for  wedges,  cones,  and  disk-shaped 
bodies. 

The  objective  of  the  present  study  is  the  numerical  investigation 
of  flows  past  inclined  plates,  especially  the  flow  characteristics  of 
vortex  shedding.  Due  to  the  high  computer  costs  involved,  the  study 
was  limited  essentially  to  the  Reynolds  numbers  10,  15,  30  and  200. 


1  References  are  listed  on  page  48= 
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FORMULATION  OF  THE  PROBLEM 


The  laminar  flow  of  an  incompressible  fluid  past  a  plate  of  elliptic 
cross-section  was  considered.  The  motion  was  assumed  to  be  two- 
dimensional  and  time-dependent.  Mathematically,  an  initial- boundary 
value  problem  for  the  Navier-Stokes  equations  must  be  solved.  This  is 
conveniently  carried  out  in  the  elliptic  coordinate  system  ( 77,8 )  (see 
Figure  1)  which  is  defined  by  the  transformation: 

2 

x  +  iy  =  a  cosh {77+  i  8),  a  >  0,  i  =  -1  .  (1) 


The  focal  distance  is  designated  by  a.  The  equations  of  motion  are 
formulated  in  terms  of  ip ,  the  dimensionless  stream  function,  and  ou, 
the  dimensionless  vorticity  component  normal  to  the  (77,  8) -plane: 


9(0 

at 


2 

Re 


2 

V  CO  , 


(2) 


(3) 


Here,  t  and  Re  are  dimensionless  time  and  Reynolds  number 

Re  =  2aU/v,  where  v  is  the  kinematic  viscosity  and  U  the  magnitude 

of  the  constant  velocity  at  infinity.  (Sometimes,  the  Reynolds  number 

Red  =  dU/v  with  d  =  2  acoshT7^  is  used  for  practical  reasons.  In 

all  cases  considered  in  this  paper,  except  for  77^  =  0.2,  Rew  Re^.) 

The  characteristic  length  and  velocity  scales  in  the  dimensionless 

quantities  are  a  and  U.  In  particular,  the  del  operator  is  made 

2  2  2 

dimensionless  by  a.  The  coefficient  h  is  defined  by  h  =cosh  77- cos  8, 
The  contour  of  the  plate  is  the  constant  coordinate  line  77=  77^ 

(Figure  1).  On  this  line  boundary  conditions  are  prescribed  such  that 
the  velocity  vector  v  is  zero.  The  dimensionless  velocity  components 


3 


const 


4 


Figure  1  -  Elliptic  coordinate  system  and  definition  of  angle  of  attack 


v  and  v.  are  connected  with  $  by 

n  d 


l  a  $ 

\  'll  38  ’ 


V  =  i  M 

e  h  a  r?  • 


Thus,  at  the  body  surface  the  boundary  conditions  are 


(4) 


n  =  7?1:  0  =  o,  1^-  =  0  .  (5) 

Here,  the  constant  value  of  0  is  chosen  to  be  zero.  Far  away  from 

-4 

the  plate,  a  uniform  parallel  flow  with  velocity  U  is  assumed,  where 
•— + 

the  vector  U  is  specified  by  a.  and  U  (Figure  1).  In  terms  of  0  the 
conditions  are: 


77 -  °° :  -  h  sin  (0  -  a),  — -  h  cos  (0-  a  )  •  (6) 

OX]  00 

Potential  flow  is  taken  to  be  the  initial  condition. 

In  order  to  compute  drag,  lift,  and  torque  on  the  plate  by  means 
of  flow  quantities  at  the  body,  the  pressure  distribution  on  the  surface 
of  the  plate  must  be  determined.  From  the  Navier-Stokes  equations  it 
follows  that 

„  -»  -2 

2  ,  -»  3  v  /V  \  ~*  /n\ 

vp  =  "  "Re  Curl  w  "  sT  "  +  v  x  00  ’  W 

where  <0  is  the  dimensionless  vorticity  vector,  and  where  the  pressure 

2 

is  made  dimensionless  by  pU  .  The  quantity  p  is  the  density  of  the 
fluid.  The  surface  pressure  is  computed  from  the  0-component  of 
(7)  to  be 


2 


5 


7F 

where  pc  is  the  pressure  at  (r/p  ^  )»  This  point  is  selected  for 

numerical  reasons.  The  pressure  pc  is  obtained  from  the  77 -component 
of  (7)  by  prescribing  the  value  of  p  at  77  =  *  to  be  p  =0: 

pc  '  J-  [  le  If  +  h  ■  hve“]  drl+3  ■  <9) 

”1  L  J 

The  drag  coefficient  is  defined  by 

Cq  -  Drag/^  a  cosh  77^  (10) 

and  consists  of  two  parts „  the  drag  coefficients  due  to  pressure  and 

friction  C_  *■  +  C__  ,  with 

D  DP  DF  ’ 

2ti  2ti 

CDp  ^ -2  tanh  t)j  cos  a  J  p^  cos  0  de  -  2  sin  a  J  p^  sine  d@  ,  (11) 

0  0 


'DF  Re 


2tt  271  : 

-  cos  01  J  co |  sin  6  d  0  +  tanh  77^  sin  a  J  60^  cos  0  d  6 


0 


0 


.  (12) 


Equation  (11)  is  simplified  by  means  of  Equation  (8)  to 

277/, 


'DP  Re 


tanhT),  cos 


2577/,  \ 

/  3co\ 
a.  t—  ) 

J  \bVL 


sin0de 


0 

277 


sin  a 


f  fr— )  COS  0  de 


Correspondingly^  the  lift  coefficient  is 
CL  ~  Lift/--  a  cosh  77^ 


(13) 


(14) 


6 


with 


4  r  27r  2tt  * 

=  ^  1  sin  a  J  tt5^sin0d6  +  tanhi7^cosa  J  u>^cos0d0  . 

e  h  o  0 


The  moment  coefficient  is 


—  2  2  2 

=  Torque/-|  U  a  cosh  t?^ 

(17) 

with 

Re  eosh\  £  S1”  ^ 

(18) 

4  277 

CMF  =  Re  tanh”l  I  “’ld6  ’ 

(19) 

0 


NUMERICAL  ANALYSIS 

The  infinite  domain  of  integration  in  the  (rj,  0) -plane  is  replaced 
by  a  finite  network  of  points  (r?^  +  (i  -  1)  At?  ,  (j  -  1) A0 )  with  i  =  1, . . .  L, 
j  =  1, . . .  M.  For  most  examples,  L  and  M  are  75  and  60;  that  is,  the 
total  number  of  grid  points  is  4500  and  A0  =  77 /30.  In  special  instances 
M  =  80  (or  A0  =  tt/40)  is  used-  For  all  cases  A 77  =  0.05.  This  choice 


7 


16  1 8 

is  based  on  previous  experience^  s  '  .  With  At?'-  0*  05  and  L  =  75s  the 
outer  boundary  n'  !  rjf  +  (L  -  1)Atj  is  about  11  plate  lengths  away 
from  the  body  center.  Figure  2  shows  the  grid  system  in  the  vicinity  of 
the  plate  for  rjj  '•  0*  1.  Two  different  layouts  are  used:  one  which 
includes  the  tip  points  in  the  grid>  another  which  excludes  the  tip  points 
by  shifting  the  grid  A 0/2  in  the  0- direction  (see  Figure  3).  The 
importance  of  this  distinction  was  noticed  in  the  course  of  the  study 
when  numerical  computations  revealed  a  much  larger  stability  for  plates 

with  high  curvature  (r?<  <  0. 1)  when  the  tip  points  were  excluded  from  the 

1  17 

grid.  A  detailed  account  was  given  in  an  earlier  report  . 

The  numerical  solution  of  the  initial  boundary  value  problem  stated 

in  the  previous  section  poses  two  essential  difficulties:  the  approximation 

of  the  differential  operators  of  Equations  (2)  and  (3)  through  suitable 

difference  schemes,,  and  the  prescription  of  boundary  conditions  for  the 

outer  boundary,,  where  <  <*  in  the  numerical  scheme* 

The  discretization  of  the  differential  operators  is  similar  to  that 
1 8 

used  by  Rimon  '  „  The  linear  part  of  Equation  (2)  is  approximated  by 

the  DuFort-Frankel  scheme*  The  Jacobian  is  written  in  conservation 

form  and  expressed  by  central  difference  formulae*  Then,  Equation  (2) 

yields,,  when  solved  for  the  (n-t-l)^  time  step  of  co,  .  , 

3 
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Figure  2  -  Elliptic  grid  system  near  the  body,  77^=0 . 1,  At)  =  0.05,  A8=ff/30.  The  total 

number  of  grid  points  is  75  x  60  =  4500 


where  the  Jacobian  is 


h2 


and  the  velocity  components 


2k  -(hV)i-i,i  j 

2ST  |(hv9“,)M+i  -  (hve“)i!j-i  j 


(21) 


(v  ) 


1  ^i,  j+1  ~  ^i,  j-1 


nj  h- 


2A6 


’  (ve\  ITT  2^  ^  •  <22> 

J  1«  .1 


It  is  recognized  that  central  difference  schemes  for  the  Jacobian  cause 
oscillations  in  oo  which  result  in  "cancer  cells".  These  oscillations 
destroy  the  quality  of  the  solutions  for  larger  Reynolds  numbers  unless 
smaller  space  increments  are  employed.  Hence,  the  capacity  of  the 
computer  sets  an  upper  limit  on  the  Reynolds  number.  With  the  space 
increments  A 77=  0.05  and  A0=  6°  generally  used  in  this  report,  the 

upper  limit  is  about  Re  =  300  (see  also  Reference  18).  Improvements 

19  11 

were  found  by  Fromm  and  Dawson  and  Marcus  who  applied  fourth 
order  approximations  and  a  combination  of  central  difference  and 
backward-forward  schemes,  respectively.  Since  the  present  study  is 
restricted  to  the  pure  K^rm^n  vortex  street  range,  with  Re  <  300,  the 
simple  central-difference  scheme  was  adopted. 

Equation  (3)  is  approximated  by  the  five-point  formula  which  yields, 
for  0.  .  , 

J 


11 


/  ...  1  _  _  JL 

Vi  i  ' '  *}  9  9 

1,1  1  (At)S  +(A9) 

-i  (At))2  .j  .+lf  .  .) -(Ar;A9)2h?  ,co,  ,1.  (23) 

Ipj  l  1?J1  l^J  IjJJ 

This  system  of  algebraic  equations  is  solved  with  Gauss-Seidel  line 
over  relaxation  applied  along  lines  of  constant  77  ,  The  overrelaxation 
factor  is  1,  82 „  The  iteration  process  is  halted  after  the  kth  iteration  if 

|vV  -  c/ |<  €  (24) 

at  each  grid  point ,  where  e  is  of  the  order  of  10'^  in  most  cases , 

For  Re  :  15  and  30,  a  45°,  e  of  the  order  10"^  had  to  be  used.  The 
number  of  iterations  varies  with  the  magnitude  of  temporal  change  of 
the  flow  field.  For  Re  "■  20Q,  «  "  45°,  77  j  “  0, 1  the  growth  of  the  vortices 
attached  behind  the  plate  requires  about  60  iterations  per  time  step, 
whereas  vortex  separation  needs  about  330, 

0 

During  the  preparation  of  this  report  several  sources  (Wirz  , 

20  21 

Briley  and  Walls  ,  Goods on*,  Buzbee  et  al  )  suggested  to  the  authors 
that  modified  Peaceman-Rachford  methods  and  methods  of  cyclic  block 
reduction  would  result  in  more  efficient  use  of  computer  time  than  the 
presented  method.  Studies  of  the  flow  problem  discussed  in  this  paper 
are  underway  to  check  the  proposed  techniques.  Results  will  be 
published  later. 


{A6)‘ 


+  *«-l  «} 


r 


*  Goodson,  R,  ,  Mechanical  Engineering  Department,  Purdue 
University,  West  Lafayette,  Indiana,  private  communication. 
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At  the  body  surface  a  one-sided  difference  scheme  must  be  used 

in  order  to  calculate  the  vorticity  a)-.  .  .  That  this  approximation 

3 

strongly  affects  the  numerical  stability  was  demonstrated  in  a  previous 
22 

investigation  „  By  trial  and  error  it  was  found  that  the  formula 


co 


1,  j 


4h?  .(At;)2 

J 


(*2,i+4*3,i'*4,i) 


(25) 


yields  the  best  results.  This  equation  is  derived  from  the  Taylor  series 

expansions  with  the  nonslip  condition  (Zty/bn)*  .  s  0  incorporated 

3 


/aiA  _  /d  ; 

'SW2,i  ‘  'i,: 


'ajs\ 

,*/l. 


(26) 


By  using  <o 


1,3 


2 

-4-  ^-4)  and  replacing  (™)  with 
h2  an2/l,i 


one  arrives  at  Equation  (25).  Hence;,  Taylor  series  expansions  are 

used  which  result  in  an  expression  for  co1  .  with  error  of  second 

order  in  Arj.  However,  when  first  derivatives  at  inner  points  are 

replaced  by  finite  differences,  the  final  formula  for  a>«  .  is  of  the 

■*■,3 

first  order.  The  hybrid  form,  Equation  (25),  is  more  stable  than  the 
first  order  approximation  used  by  most  other  investigators 
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+  O  (A  ri ) 


(28) 


h2  .(At?)2 
A5  J 


and  much  more  stable  than  second  order  approximations  investigated, 
for  instance. 


2hj 

■*'$  J 


(29) 


In  fact,  it  was  found  that  At  is  about  10  times  smaller  for 
s  max 

Equation  (28)  and  100  times  smaller  for  Equation  (29)  than  Atmax  for 

Equation  (25) „  However,  the  statement  that  Equation  (25)  is  superior 

23 

to  Equation  (28)  is  not  valid  in  general,,  For  instance,  in  other  work 
we  found  Equations  (25)  and  (28)  to  be  comparable  with  regard  to 
stability,, 

Since  the  domain  of  numerical  integration  is  bounded,  severe 
difficulties  are  encountered  in  prescribing  the  conditions  at  the  outer 
boundary  77"  r?L  .  This  boundary  is  almost  circular  and  is  divided 
in  half,  On  the  upstream  half  of  the  boundary,  the  conditions  (6)  are 
used  with  the  alteration  that  the  second  condition  is  replaced  by 
vanishing  vorticity: 

0  £  6  £  |  +  a 

V “ Vt  ,  3  :  M,  h  sin(e "  <*)>  «=  0  •  (30) 

L  |ir+a  ^  6  <  2tt  drj 

In  difference  form  (tip  point  not  included  in  grid): 

fekrV “l,s=0-  (31> 
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On  the  downstream  half  of  the  boundary  a  convection  type  condition  is 
used: 

7T  3 

77=  r?L  ?  2  +  oc<'Q<'21I  +  oi  : 


5  co 
at 


+ 


4(U*  v)co  =  0  , 


at  +  u 


av 


(32) 


This  type  of  condition  was  proposed  by  Dawson  and  Marcus  for  the 
vorticity  and  seems  to  work  well  when  extended  as  above  to  the  velocity. 
With  the  difference  approximations  of  Equations  (32)  the  values  of  co 
and  3^/377  are  computed  on  the  boundary  77=  7?T  : 


n+1  n  At 
wL,rwL,J  +  ll2 


L,j 


(cos  a  sinh  ?7t  cos  0. 


(n  n 
^L,  j  "  wL-15  j\ 

- U - Li  ! 

At?  I 


+(sina  sinhr?T  cos  0.  -  cos  a  cosh77T  sin  0.) 

-Lj  J  J 
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(M\n+1 ,  Mn  +  At 
Wt,  i  V 


^r)/Lj  WlJ  2^e 


(cos  a  sinhr?L  cos  0^  +  sin  a  cosh  77^  sin  @j+j.) 


n 


”L  ,j+l 
hL,  j+1 


-  (cos  a  sinhr^cos  0.  ^  +  sin  a  coshr)^ 


n 


V 


sin  0. 


L,j-1 


H  Vh 


v 


4  (sin  a  sinhr?^eos  0  ^  -  cos  acoshr?^  sin  0.+ j)  — ■■ 


©Lj  j+1 


J+i  “L,  j+l. 


(sino;  sinh?7L  cos  0.  ^  -  cos  accsh77L  sin  0._.j)  ^ 


e,  j-i 
L3  j- 1 


+  2A0  (cos  a  sinhr?T  cos  0.  +  sinacoshr7T  sin  0.)oOj  . 

ij  ]  Cj  J  Ijj  j 


(34) 


The  decision  to  use  convection  equations  of  the  type  in  Equations  (32) 
is  heuristic  and  guided  by  the  idea  of  interfering  with  the  downstream 
flow  as  little  as  possible.  Photographs  of  wind-tunnel  flows  reveal  that 
the  shapes  of  the  vortices  are  well  preserved  over  many  cycles  of  the 
Karman  vortex  street. 

The  integration  process  is  carried  out  in  the  following  way:  The 

yi  |  1 

vorticity  00  -  for  the  advanced  time  step  (n+1)  is  computed  at  the 

inner  points  according  to  Equation  (20).  Next,,  o>,  .  and  ( S ^  / 5 77  )T  . 

J  1-Jj  J 

are  determined  on  the  downstream  boundary  points  by  means  of 
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Equations  (33)  and  (34).  Then  .  is  calculated  with  the  aid  of 
Equation  (23).  The  cycle  concludes  with  the  calculation  of  . 
from  Equation  (25). 

The  maximum  stable  time  step  At  ,  beyond  which  numerical 
instability  occurs ,  is  determined  by  increasing  the  time  step  until 
oscillations  from  one  time  step  to  another  appear  in  the  co-values 
near  the  tips.  The  magnitude  of  Atmgx  depends  on  factors  other 
than  co^  .  It  is  a  linear  function  of  the  Reynolds  number,  at  least  for 

Re  £  200,  and  decreases  rapidly  as  either  A?),  A0,  or  ry1  approaches 

18  ^ 
zero  .  There  seems  to  be  no  dependence  on  a.  Computer  time  is 

saved  by  enlarging  77^  as  much  as  possible  without  the  loss  of  the 

essential  features  of  the  flat  plate  (if  77^  is  changed  from  0.05  to  0. 1, 

At  increases  fourfold).  The  values  of  At  are  for  771  =  0. 1,  tip 
max  max  l 

point  out, 


Re 

a 

At 

max 

15 

Ac0 

45 

0.0004 

30 

45° 

0.0007 

50 

90° 

0.0012 

200 

0°,  45° 

0.005 

Smaller  time  steps  are  necessary  near  t  =  0,  when  large  vorticity 

gradients  are  present.  For  high  surface  curvature,  that  is  for 

77^  <  0. 1,  stability  is  improved  if  the  grid  system  is  shifted  by  A0/2 

in  the  0-direction  so  that  there  are  no  grid  points  at  the  tips.  This 

does  not  affect  the  overall  accuracy  of  the  solution  and,  in  addition  to 

improving  the  stability,  makes  possible  the  computation  of  the  flow 
17 

field  when  77^  =  0  . 
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Accuracy  considerations  are  similar  to  those  for  the  flow  past 
18 

oblate  spheroids  „  Additional  information  on  the  quality  of  the  solution 

is  obtained  by  integrating  do o/brj  over  the  body  contour,,  The  integral 

of  this  function  over  a  segment  of  the  surface  is  the  pressure 

difference  between  the  end  points  of  the  segment  Hence^  the 

integration  around  the  entire  contour  of  the  body  should  yield  zero 

and  is  a  sensitive  test  which  must  be  passed  to  a  certain  degree  by  any 

o 

accurate  solution,,  For  Re  20Q„  a  ?-■  45  ,  "  0, 1,  the  value  of  this 

integral  is  always  smaller  than  3%  of  the  difference  between  the 
maximum  and  minimum  pressure  at  the  surface,,  Other  ’'confidence- 
building''  tests  compare  our  results  to  those  obtained  with  different 
methods . 

The  numerical  calculation  of  the  drag?  lift„  and  moment  coefficients 

is  sensitive  to  the  difference  scheme  used  for  the  approximation  of 

1 7 

the  integrals  of  Equations  pi)  through  (19)  »  These  coefficients  can 

be  computed  either  with  the  aid  of  a  double  integral  (for  C^p  it  is 
formula  (11)  )  or  by  means  of  a  single  integration  (for  Cpp  it  is 

formula  (13)  )„  For  both  methods  two  difference  schemes  have  been 

.  2tz 

checked  and  are  presented  here  for  the  example  p.  cos  0d  0  with 

0  1 

the  tip  not  included  in  the  grid^ 


Simpson's  rule” 


J  p  j  cos  0  d  0 
0 


A0 

3 


Aft  3 

2  p  t  j  cos  4  4  P|  2  cos  ^  2  A  0  ^ 


5  2M-1 

2  Pj  g  cos  (g  A0)  t...+4pj  Mcos(--~2“— A0) 


(35) 
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Second  approximation: 

2^  M  /  1  i 

J*  P1cos0d0=  E  A0  P1  jCOS^(j-g)A0  .  (36) 

0  3  ^ 

For  the  grid  which  does  not  contain  the  tip  point  the  following  values 
have  been  computed:  Re  =  50,  a=  90°,  r}^  =  0.1,  t  =  1.71. 


CDP 

CDF 

Double  integration: 

Simpson's  rule: 

5.67 

0.230 

Approximation  (36): 

5.67 

0.230 

Single  integration: 

Simpson's  rule: 

6.50 

0.230 

Approximation  (36): 

5.67 

0.230 

If  the  tip  point  is  included  in  the  grid,  the  data  for  the  case  Re  =  50, 
a-  90°,  r /j  -  0.05,  t  =  2.31  are 


Double  integration: 


Simpson's  rule: 

5.21 

0.097 

Single  integration: 

Simpson's  rule: 

2.70 

0.097 

Approximation  (36): 

5.20 

0.105 

17 

Based  on  earlier  results  where  theory  and  experiment  were  compared, 
it  is  assumed  that  the  worst  result  is  obtained  with  the  single  integration 
by  using  Simpson's  rule.  Unfortunately,  this  case  was  reported  in  the 
condensed  version  of  this  paper,  Lecture  Notes  in  Physics  8,  1971, 
Springer  Verlag,  p.  78.  The  value  of  2.  8  given  in  this  paper  for  Re  =  50, 
a=  90°,  =  0.05,  t  =  2..13  is  not  correct.  Throughout  the  present 

report  the  single  integration  by  using  approximation  (36)  is  employed. 
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RESULTS 


All  computations  were  carried  out  in  double  precision  on  an 
IBM  360-91  computer,.  The  graphic  display  of  streamlines  and  lines 
of  constant  vorticity  was  produced  with  a  Stromberg-Carlson  SC-4020 
charactron  plotter,,  The  examples  selected  for  computation  are 
compiled  in  Table  1. 

Symmetric  Flows 

For  the  symmetric  flow  configurations ,  which  according  to 
experiments  are  expected  for  a."  0°,  Re  <  200,  77^  s  0„  1  and  a.  =  90°, 

Re  <  50?  77 ^  arbitrary,  a  steady  state  is  approached  as  t  os  .  For 
Re  -  10,  a  90  ,  77^  “  0.  1  one  aspect  of  the  transient  stage  is  displayed 
in  Figure  4.  where  the  drag  coefficient  is  plotted  against  time.  Figures  5 
and  6  show  the  streamlines  and  lines  of  constant  vorticity  for  the  almost 
steady  state  t  -  11,2  ,  The  development  of  the  twin  vortices  behind  the 
body  is  recorded  in  Figure  7,  where  the  dimensionless  length  of  the 

wake  s/d  is  plotted  against  time.  Although  the  steady-state  value 

24 

measured  by  Taneda  is  approached  quite  accurately,  the  data  for  the 

transient  stage  are  low  compared  with  the  empirical  formula  developed 

25 

by  Taneda  and  Honji 

s/d  -  0o  89  (t/2)2/3,  Re  *  18. 1  (37) 

This  relation,  which  does  not  contain  the  kinematic  viscosity  v , 

holds  only  for  an  intermediate  time  interval.  As  can  be  seen  immediately 

from  Figure  7,  the  wake  length  approaches  a  finite  value  when  t  -*  » , 

a  process  which  is  governed  by  viscous  forces.  Near  t  =  0,  Equation  (37) 

25 

cannot  be  valid  either.  As  shown  also  in  Taneda  and  Honji’s  experiments  , 
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TABLE  1  -  COMPILATION  OF  THE  CALCULATED  EXAMPLES 


Re  , 
d 

a 

*1 

1  FINAL 

GRID 

TIP  IS  GRID  POINT 

10 

90° 

0. 1 

11.2 

75  x  60 

NO 

15 

45° 

0o  1 

7.2 

75  x  60 

NO 

30 

45° 

0.1 

16.3 

75  x  60 

NO 

50 

0° 

0.  1 

1.5 

75  x  60 

YES 

50 

to 

o 

o 

0.05 

2.3 

75  x  60 

YES 

50 

90° 

0.1 

3.5 

75  x  60 

NO 

200 

0° 

0. 1 

18.7 

75  x  60 

YES 

200 

o 

o 

0 

1.9 

75  x  80 

NO 

200 

45° 

0 

0.04 

75  x  80 

NO 

200 

45° 

0.1 

30.0 

75  x  60 

NO 

200 

45° 

0.2 

8.1 

75  x  60 

NO 

200 

4* 

CJl 

o 

0.2 

11.5 

75  x  80 

NO 

21 
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Figure  4  -  Drag  coefficient  versus  time  for  Re 


Figure  5  -  Streamlines  for  Re  =  10,  a  *  90°,  tn  =  0. 1  at  t  *  11.2 

(almost  steady  state) 
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Figure  6  -  Lines  of  constant  vorticity  for  Re  =  10,  a*  90°, 
772  =  0.1  at  t  =11.2  (almost  steady  state) 
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the  vortices  are  generated  at  the  tips  and  then  join  at  the  centerline 
after  a  certain  time  (Figure  8)  .  The  instant  at  which  they  join  is  the 
beginning  of  the  wake  development  with  nonvanishing  s.  For  small 
Reynolds  numbers  the  two  time  periods  which  are  excluded  from 
Equation  (37)  overlap  and  the  curve  of  Equation  (37)  is  not  reached. 
This  appears  to  be  the  reason  for  the  low  values  at  Re  ~  10„ 

The  drag  coefficients  for  the  almost  steady  state  are  compared 
with  other  sources  in  Table  2. 


TABLE  2  -  DRAG  COEFFICIENTS  FOR  SYMMETRIC 
FLOW  CONFIGURATIONS 


O' 

Re 

*1 

1  FINAL 

CD 

Computed 

4 

From  Dennis  &  Chang 

■■ 

50 

0. 1 

1.50 

1.39 

1.432  (for  Re  =  40) 

200 

0.1 

18.7 

0.55 

0.  492  (for  rjj  -  0) 

m 

10 

0. 1 

11.2 

6.02 

The  values  for  a  ~  0°  compare  favorably  with  the  steady-state  solutions 

4 

of  Dennis  and  Chang  . 

It  may  be  mentioned  that  in  the  condensed  version  of  this  paper, 

published  in  the  Lecture  Notes  in  Physics,  No.  8,  Springer -Verlag,  1971, 

the  case  a--  90°,  Re  =  50,  r\ j  -  0.05,  t  ~  2. 13  was  considered  as  an 

almost  steady  state.  This  statement  is  incorrect,  since  the  wake 

24 

grows  considerably  larger  when  t  „ 
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Figure  8  -  Sequence  of  streamlines  in  the  initial  phase  of  wake 
development  for  Re  =  50,  a=90°,  =0. 1 


Nonsymmetric  Flows 

Most  of  the  numerical  data  compiled  are  obtained  for  the  non¬ 
symmetric  case  of"  45°,  Re  =  200.  The  first  series  of  pictures,  which 
is  shown  in  Figure  9,  is  selected  to  display  the  flow  behavior  following 
the  abrupt  start  of  the  body  at  t  =  0.  Immediately  after  t  =  0,  in  the 
initial  potential-flow  field,  a  starting  vortex  forms  and  separates  from 
the  trailing  edge.  In  the  streamline  pictures  this  vortex  is  visible  as 
a  wavy  pattern  which  indicates  the  movement  of  the  vortex  relative  to  the 
plate  (which  is  kept  at  rest).  The  lines  of  constant  vorticity  clearly 
exhibit  a  local  extremum  at  the  center  of  the  vortex.  In  the  early  time 
period  the  zero  streamline  migrates  from  the  rear  stagnation  point  of 
the  potential  flow  to  the  trailing  edge.  This  verifies,  even  for  the  low 
Reynolds  number  Re  -  200,  the  hypothesis  of  Lanchester-Prandtl,  that 
circulation  around  an  airfoil  is  generated  by  a  starting  vortex  and  that  the 
zero  streamline  follows  a  kind  of  Kutta  condition. 

The  computation  of  the  case  Re  ~  200,  o:  =  45°,  ~  0. 1  was 

carried  out  over  three  cycles  of  vortex  shedding.  The  last  cycle  was 
repeated  with  the  fine  mesh  A rj~  0.05,  A8  =  7t/40.  The  flow  patterns  of 
the  third  cycle  obtained  with  the  fine  mesh  are  displayed  in  Figure  10. 

The  development  and  the  detachment  of  the  vortices  are  similar  to  those 
of  the  first  cycle  in  Figure  9.  An  overall  view  of  the  wake  is  presented 
in  Figure  11.  To  the  right  the  boundary  of  the  grid  system  is  visible. 

The  location  of  the  vortices  shows  that  the  vortex  street  is  not  parallel 
to  the  undisturbed  flow. 

In  Figure  12  the  drag,  lift,  and  torque  coefficients  are  plotted 
against  time  t.  The  abrupt  start  of  the  body  requires  infinite  values 
for  CD  and  C^,  whereas  the  C^-value  begins  from  zero.  The  Strouhal 
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Figure  10  -  Sequence  of  streamlines  and  equal-vorticity  lines  for  the 
third  cycle  of  Re  =  200,  a=  45°,  r]-,  =  0. 1  after  the  abrupt  start 
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Figure  12  -  Drag,  lift,  and  moment  coefficients  versus  time  for 


number  which  is  defined  by 


(38) 


is  about  0.  23  for  the  second  cycle  and  0.25  for  the  third.  In 
Equation  (38)  n  is  the  frequency  of  the  vortex  shedding.  If  one 
relates  the  Strouhal  number  to  the  projected  plate  width  d  sin  a,  the 
values  for  the  second  and  third  cycles  are  0. 163  and  0. 177,  respectively. 

A  comparison  of  Figures  10  and  12  reveals  that  the  coefficients 
Cqj  Cl,  and  assume  relative  maxima  whenever  a  vortex  separates 
from  the  leading  edge. 

In  Figure  13  the  coefficients  CDF  and  C^F  are  plotted  against 
time  for  Re  =  200,  a  =  45°,  r\ ^  =  0. 1. 

For  =  0.2  computations  have  been  made  with  A0  =  tt/30  as  well 
as  with  A0  =  tt/40.  The  differences  in  the  CQ,  CL,  and  values 
between  the  two  cases  are  so  small  that  they  cannot  be  observed  on  a 
graphic  display. 

In  Figures  14  and  15  the  surface  vorticity  and  the  surface  pressure 
are  plotted  against  0  for  the  third  cycle  of  the  case  Re  =  200,  cl-  45°, 
r)^  =  0. 1  computed  with  the  fine  mesh.  The  production  of  vorticity  at 
the  edges  is  smallest  when  a  vortex  separates  from  the  edge,  and  it 

is  highest  when  a  vortex  starts  growing. 

26 

Based  on  Timme's  experiments  the  decay  of  vortices  in  the 
Karman-vortex  street  can  be  described  by  the  Hamel-Oseen  solution 
for  each  vortex: 
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Figure  13  -  Coefficients  for  frictional  drag  and  lift  versus  time  for 
Re  =  200,  a=  45°,  Vl  =  0. 1 


Figure  14  -  Surface  vorticity  versus  0  for  the  third  cycle  of 
Re  =  200,  a  =  45°,  ^  =  0. 1 


Figure  15  -  Surface  pressure  versus  0  for  the  third  cycle  of 
Re  =  200,  a=  45°,  r/1  =  0. 1 


where  the  polar  coordinates  (r,  <j> )  with  the  corresponding  velocity 

components  (v^  ,  v^)  are  used.  Equation  (39)  represents  the  decay 

of  a  potential  vortex  from  the  time  t*  -  t'  .  This  expression 

0  27 

represents  the  lowest  mode  of  plane  disturbances .  The  vorticity  at 
the  vortex  center  (where  the  extremum  of  vorticity  is  located)  dissipates 
according  to 


,  const 

Wr-0  "  2i/(t*-tj) 


(40) 


For  the  next  higher  mode  of  the  spectrum  of  disturbances  the  vorticity 
would  decay  as  (t'-t^)  ,  as  was  shown  elsewhere  „  The 

numerical  data  can  be  checked  against  the  analytic  prediction.  In 
Figure  16  the  decay  of  the  central  vorticity  of  the  initial  vortex  for 
Re  =  200,  a~  45  ,  77 j  -  0. 1  is  displayed  with  logarithmic  scales.  The 
best  fit  to  a  linear  relation  between  the  logarithms  of  the  central 
vorticity  and  time  is  obtained  when  t^  -  -0.3  and  the  slope  is  -1. 

Vortices  which  are  generated  after  the  initial  vortex  behave  similarly. 
Near  the  outer  boundary  of  the  grid,  where  the  mesh  size  increases 
rapidly ,  the  numerical  values  deviate  from  the  analytic  curve.  This 
deviation  depends  on  the  mesh  size  (see  Figure  16).  Apparently,  the 
numerical  solution  is  not  accurate  enough  to  describe  the  vortex  decay 
in  that  part  of  the  grid. 

o 

For  Re  =  200,  some  information  was  obtained  when  a~  10  ,  rj^-O. 
The  initial  phase  up  to  t  =  1.  94  is  recorded  in  Figures  17  and  18.  Non¬ 
periodic  vortex  shedding  from  the  leading  edge,  which  has  been  observed* 


*  Laitone,  E.V.,  University  of  California,  Berkeley,  private 
communication. 
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Figure  16  -  Decay  of  central  vorticity  of  the  initial  vortex  for 
Re  *  200,  a  *=  45°,  *=  0. 1 

-  I 
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40 


for  this  small  a -value  but  at  much  higher  Reynolds  number,  does 
not  occur  in  the  present  example,, 

For  Re  =  30,  a  =  45°,  -  0.1  numerical  results  are  obtained 

up  to  t  =  16.  Neither  a  clear  development  of  a  Karman  vortex  street 
nor  a  tendency  towards  a  steady  state  can  be  observed.  The  streamline 
patterns  in  Figure  19  show  vortex  shedding.  However,  the  corresponding 
lines  of  constant  vorticity  do  not  exhibit  local  extrema;  they  only 
oscillate.  Apparently,  diffusion  of  vorticity  at  this  low  Reynolds 
number  is  already  so  dominant  that  vorticity  extrema  do  not  develop. 
Whether  an  oscillating  state  or  a  steady  state  is  approached  when  t  ->  <» 
cannot  be  determined.  However,  the  critical  Reynolds  number  which 
separates  these  two  states  must  be  close  to  Re  =  30.  (For  bodies 
normal  to  the  flow  the  critical  Reynolds  number  is  higher  and  is  about 
45.)  Figure  20  shows  the  coefficients  CDF  and  CLF  as  functions  of 
time. 

The  case  Re  =  15,  a.  -  45°,  r?-^  -  0. 1  is  computed  up  to  t  =  7. 

Although  a  steady  state  is  expected  for  t  ,  the  time  t  =  7  is  too 
small  to  be  considered  as  giving  an  almost-steady  state.  At  least  a 
qualitative  picture  of  the  flow  patterns  is  displayed  in  Figure  21. 
Additional  information  is  recorded  in  Figures  22  and  23. 


CONCLUSIONS 

Computer  simulation  of  vortex  and  vorticity  shedding  is  possible 
with  accuracy  of  the  order  of  water-  or  wind-tunnel  experiments  for 
two-dimensional  flows  around  bodies.  With  the  space  increments  and 
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Figure  20  -  Coeffici 


Lents  versus  time  for  Re  =  15 
abrupt  start 
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number  of  grid  points  (75  x  60  or  75  x  80)  used  in  this  paper,  the 
Reynolds  number  is  limited  to  about  300.  The  flow  region  far  from 
the  body  is  not  well  represented.  Improvements  must  be  made  in 
accuracy  and  in  economy  of  computer  time.  Promising  studies  are 
under  way. 

The  sudden  start  of  a  thin  elliptic  cylinder  under  an  oblique  angle 
of  attack  produces  an  initial  vortex  which  separates  from  the  trailing 
edge.  After  this  time  the  zero  streamline  of  the  wake  leaves  the  body 
parallel  to  the  major  axis  at  the  trailing  edge.  This  phenomenon, 
although  occurring  at  Re  =  200,  verifies  the  soundness  of  the  Lanchester- 
Prandtl  hypothesis  of  potential  flow,  in  which  a  circulation,  necessary 
for  the  existence  of  lift,  is  postulated  and  determined  by  means  of  the 
Kutta  condition. 

For  a  thin  elliptic  cylinder  with  a  -  45°  the  lower  limit  of  vortex 
shedding  appears  to  be  about  Re  =  30.  (For  symmetric  configurations 
normal  to  the  flow  the  lower  limit  is  about  Re  =  45. ) 

All  vortices  shed  from  the  body  dissipate  in  such  a  way  that  the 
relative  extrema  of  their  vorticity  (which  may  be  considered  as  vortex 
centers)  decrease  in  value  approximately  as  l/(t- const).  Hence, 
they  behave  like  decaying  potential  vortices  when  separated  from  the  body. 
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